Whole‐mitogenome analysis unveils previously undescribed genetic diversity in cane toads across their invasion trajectory

Abstract Invasive species offer insights into rapid adaptation to novel environments. The iconic cane toad (Rhinella marina) is an excellent model for studying rapid adaptation during invasion. Previous research using the mitochondrial NADH dehydrogenase 3 (ND3) gene in Hawai'ian and Australian invasive populations found a single haplotype, indicating an extreme genetic bottleneck following introduction. Nuclear genetic diversity also exhibited reductions across the genome in these two populations. Here, we investigated the mitochondrial genomics of cane toads across this invasion trajectory. We created the first reference mitochondrial genome for this species using long‐read sequence data. We combined whole‐genome resequencing data of 15 toads with published transcriptomic data of 125 individuals to construct nearly complete mitochondrial genomes from the native (French Guiana) and introduced (Hawai'i and Australia) ranges for population genomic analyses. In agreement with previous investigations of these populations, we identified genetic bottlenecks in both Hawai'ian and Australian introduced populations, alongside evidence of population expansion in the invasive ranges. Although mitochondrial genetic diversity in introduced populations was reduced, our results revealed that it had been underestimated: we identified 45 mitochondrial haplotypes in Hawai'ian and Australian samples, none of which were found in the native range. Additionally, we identified two distinct groups of haplotypes from the native range, separated by a minimum of 110 base pairs (0.6%). These findings enhance our understanding of how invasion has shaped the genetic landscape of this species.


| INTRODUC TI ON
Invasive species are a substantial global concern, particularly due to likely future range expansions under climate change (Bonebrake et al., 2018;Seebens et al., 2017), an increasing rate of anthropogenic introductions (Lowe et al., 2000), and their negative impacts on native biodiversity and the economy (Bradshaw et al., 2016;Mainka & Howard, 2010).Biological invasions sometimes present a genetic paradox where, despite small founding population sizes that often reduce genetic diversity in the introduced range (i.e.genetic bottlenecks), they still maintain the ability to readily adapt to novel conditions (Allendorf & Lundquist, 2003;Estoup et al., 2016;Schrieber & Lachmuth, 2017).Many invasive populations have shown evidence of adaptation to their new environments (Rollins et al., 2013), creating an opportunity to study the process of rapid evolution.
Various genetic markers, including those from the nuclear and mitochondrial genomes, have been employed to investigate the introduction history and evolutionary processes of invasive species (McGaughran et al., 2024).Because of their distinct mode of inheritance and mutation rates, nuclear and mitochondrial markers may yield different insights into the evolutionary events underlying invasion (Toews & Brelsford, 2012).Genetic variation in invasive species is shaped by complex interplays among genetic drift, demography, and selection.Small founding populations, common during an invasion, can promote genetic drift, leading to reduced genetic diversity.Furthermore, while genetic drift results in less efficient selection (Gravel, 2016), populations expanding their range into new environments may experience novel selection regimes (Sakai et al., 2001) and spatial sorting (Shine et al., 2011) that further reduce genetic diversity.
Mitochondrial DNA (mtDNA) is more susceptible to genetic drift than nuclear DNA due to the former's uniparental inheritance and smaller effective population size (Klucnika & Ma, 2019;Xiao et al., 2017).Nonetheless, mtDNA has been extensively used to study the evolutionary history of different species or populations within the same species because of its non-recombining nature and purported near-neutrality.However, the neutrality of mtDNA has been challenged by evidence of selection, likely due to the important role of mitochondria in ATP production and cellular activities (Meiklejohn et al., 2007).Deleterious mutations are rapidly removed by purifying selection to maintain the function of the electron transport chain.Furthermore, variation in mtDNA has been linked to adaptation to high latitudes (Chen et al., 2020;Wang et al., 2021) and temperature (Baker et al., 2019;Cheng et al., 2013;Silva et al., 2014).
Despite these challenges, there exists a substantial repository of mtDNA data in public databases, enabling comprehensive comparative genetic analyses, increasing its usefulness as a marker.
The cane toad (Rhinella marina) is one of the world's most infamous, and well-studied, invasive species (Lowe et al., 2000).Cane toads are native to Central and South America (Zug & Zug, 1979) and have been introduced to multiple Caribbean islands for biological control of cane beetles, a parasite of sugarcane (Easteal, 1981).
Cane toads were transported from the Caribbean to Hawai'i in 1932 for cane beetle control, and in 1935, 101 Hawai'ian cane toads were translocated from Oahu to Queensland, Australia (Turvey, 2013).
However, they failed to control cane beetles in Australia and, facilitated by a high reproductive rate, subsequently spread across the continent to arid Western Australia.This species' inexorable spread throughout its introduced range in Australia has had catastrophic effects on the naïve predators it has encountered (Chinchio et al., 2020).Freshwater crocodiles, snakes, lizards, and quolls have been particularly badly affected by the cane toad invasion in Australia (Shine, 2010).The colonisation of the cane toad is still ongoing; toads occupy over 1.2 million km 2 of Australia (Urban et al., 2008) including tropical Queensland, the Northern Territory, and the Kimberley region of northern Western Australia (Radford et al., 2019).
Like in many other invasions, the cane toads' translocation resulted in genetic bottlenecks.Australian cane toads have very low major histocompatibility complex (MHC) class I and class II (Lillie et al., 2014(Lillie et al., , 2016(Lillie et al., , 2017) ) and microsatellite (Leblois et al., 2000) diversity.Previous analysis of the NADH dehydrogenase 3 (ND3) mitochondrial gene from 31 individuals sampled from Hawai'i and Australia identified a single haplotype, a substantial reduction of haplotype diversity compared to the five haplotypes observed in the native range (Slade & Moritz, 1998), further supporting the presence of genetic bottlenecks in the invasive ranges.Despite this low genetic variation, toads from the Australian range differ substantially in immune function, dispersal ability, and behaviour (Alford et al., 2009;Brown, Phillips, Dubey, & Shine, 2015;Gruber et al., 2018;Gruber, Whiting, Brown, & Shine, 2017;Phillips et al., 2010).Some of these traits have been shown to be intergenerationally transmitted, consistent with genetic causes.For example, toads from range-core and range-edge populations that were raised in a common-garden experiment showed divergence of behavioural traits and morphology (Gruber, Brown, Whiting, & Shine, 2017;Hudson, Brown, & Shine, 2016).The invasion success of this species in Australia, coupled with evidence of evolutionary change since introduction, makes it an ideal system for the study of rapid evolution during invasion.
Despite extensive research into toad evolution in the introduced range, we need improved genetic resources to clarify the molecular mechanisms underlying these changes.Transcriptome data from spleen, muscle, and/or brain tissues of toads from French Guiana, Hawai'i, and Australia and liver tissue from Brazil exist (Richardson et al., 2018;Rollins et al., 2015;Selechnik, Richardson, Shine, Brown, & Rollins, 2019;Yagound et al., 2022), and a draft reference genome (Edwards et al., 2018) is available, but population genomic analyses of cane toads across this invasion trajectory have thus far been restricted to a single study of single-nucleotide polymorphisms (SNPs) from the nuclear genome (Selechnik, Richardson, Shine, DeVore, et al., 2019).Here, we present the first complete cane toad mitochondrial genome (mitogenome), assembled from the same Western Australian individual as the draft nuclear genome (Edwards et al., 2018).We then investigate evolution across the invasion trajectory using mtDNA data derived from RNA-Seq and wholegenome resequencing (WGS).We compare a subset of our data to those of a previously published study of ND3 on the same populations (Slade & Moritz, 1998), to determine whether our wholemitochondrial genome analysis can extend population genetic inferences.We predicted that there would be decreasing genetic diversity within the mitochondrial genome with increasing distance from the native range due to genetic bottlenecks.We also predicted that we would identify more genetic diversity within Australia than has been previously characterised, given the phenotypic evidence of adaptation to this novel environment.

| Reference mitochondrial genome
High-molecular-weight genomic DNA was extracted from the liver of a female cane toad collected in the Kimberley region of Western Australia and sequenced on the Illumina HiSeq X Ten and PacBio RS II platforms as described previously (Edwards et al., 2018).We assembled the mitogenome from 36 SMRT cells generated from the first two SMRTBell libraries from Edwards et al. (2018) (ENA experiments ERX2389178 and ERX2389179).
We downloaded the 17,757 base pair (bp) Bufo japonicus complete mitochondrial genome (gi: 157939553) from NCBI on 29 March 2016 and used it to identify mitochondrial subreads from the raw sequencing data.We used GABLAM (Davey et al., 2006) to map the B. japonicus mtDNA onto cane toad PacBio subreads using blastn (blast+ v2.2.31) and extracted 143 subreads that had 75% + bidirectional coverage (e.g. at least 75% of the subread mapped onto the mtDNA and at least 75% of the mtDNA was covered by the subread).We reverse complemented any subreads that mapped the negative strand.We identified the subread with the best global hit to the B. japonicus mtDNA and performed a second GABLAM search of all 143 subreads against this new query, retaining 29 sequences that had >80% global sequence identity within local matches in the same orientation and order, excluding any "wrapping" due to circularity.We aligned these 29 sequences using MAFFT v2.273 (Katoh & Standley, 2013) with a very small gap penalty to account for the high indel rate (--localpair --maxiterate 1000 --op 0.1 --ep 0.1 --lop -0.1 --nuc).Using SeqSuite (Edwards et al., 2020), we generated a consensus sequence from the alignment using the most frequent non-gap base call at each position.We excluded from the consensus sequence any column with fewer than 5 non-gap sequences.We mapped all 143 raw subreads onto the consensus sequence using BLASR and polished it using Quiver (SMRT Analysis v2.3.0).We used GABLAM to identify the overlapping ends of the sequence, which was then circularised by trimming off the first 1957 bp for a final length of 18,152 bp.We performed a final polishing step using the Illumina data (ENA experiment ERX2845325), mapped with Bowtie v2.3.3 (Langmead & Salzberg, 2012) and error-corrected with Pilon v1.20 (Walker et al., 2014).
We annotated the mitochondrial genome using MITOS online (Bernt et al., 2013) with the vertebrate mtDNA genetic code (Genetic Code Table 02).We manually curated the annotation of protein-coding genes and associated stop codons where required: where no canonical in-frame stop codon was found, we assumed that a partial (T or TA) stop codon would be completed upon polyadenylation.We compared the gene arrangement of the mitochondrial genome to the closely related Rhinella species (Rhinella cf.acrolopha) (accession number: KT221613) (Machado et al., 2016).We calculated codon usage using https:// www.bioin forma tics.org/ sms2/ codon_ usage.html (accessed on 29 March 2023).We uploaded the final annotated mitochondrial genome to NCBI (accession number: NC_066225.1).We drew the mitogenome map using OGDRAW (Greiner et al., 2019).

| Sequence data for population-and species-level analyses
We downloaded previously published cane toad RNA-Seq data from different tissues from NCBI BioProject PRJNA510261 and PRJNA395127 (spleen), PRJNA479937 (brain), and PRJNA277985 (muscle) for analysis.We included a total of 125 individuals comprising 144 tissue samples, from individuals sampled in French Guiana (n = 24), Hawai'i (n = 8), and Australia (n = 93) (Figure 1; Table S1).Where multiple tissues were available for an individual (brain and spleen), we used the more common tissue type (brain).
We produced WGS data for an additional 15 cane toad samples (see "WGS" in Table S1) to provide a broader geographical context in the native range and increase our whole-mitochondrial genome population genetic analysis sample size, bringing our total number of individuals included to 132.Individuals were chosen where possible to enable comparison of DNA-Seq-and RNA-Seq-derived mitogenomes (8 individuals).We extracted genomic DNA using a Gentra PureGene DNA extraction kit (QIAGEN) (Figure 1, Table S1).For two of the samples (RMH006 and RM265), paired-end sequencing libraries were prepared using the TruSeq DNA PCR-Free kit (Illumina) and then sequenced on the Illumina NovaSeq 6000 platform at the Ramaciotti Centre for Genomics, University of New South Wales, Sydney, Australia (UNSW Sydney).
For the other 13 samples, paired-end sequencing libraries were prepared according to the manufacturer's specifications for the NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs).Three of these (samples RMF031, RMF042, and RMF044) were then sequenced on the 10X Genomics Chromium platform at the Ramaciotti Centre for Genomics, UNSW Sydney.The remaining 10 libraries were sequenced on the Illumina NovaSeq 6000 platform at Deakin Genomics Centre, Victoria, Australia, with a target depth of 20-fold coverage per genome.We processed raw sequencing reads as described above.
To put our data into the context of previous studies (Acevedo et al., 2016;Slade & Moritz, 1998), we obtained the mitochondrial ND3 gene with tRNA-Gly and tRNA-Arg flanking sequence (ND3 dataset hereafter) data from NCBI (Table S1).We aligned our data with this ND3 dataset using the R package DECIPHER (Wright, 2016) and trimmed our data to match the pre-existing sequences.

| Construction of mtDNA genome and variant calling
We used GATK HaplotypeCaller (Poplin et al., 2018) to call SNPs and indels using the following criteria: soft clipped bases were not included, and variants with minimum Phred-scaled confidence of 20 threshold were kept.To filter low-quality variants, we used GATK VariantFiltration (Van der Auwera, 2020) to remove variants in the resulting VCF file: clusters of 3 SNPs within a window of 35 bases, variants with FisherStrand (FS) >30.0,QualByDepth (QD) <2.0, depth of coverage (DP) <20.0, and allelic frequency (AF) <0.05.Due to the difficulty of mapping short-read RNA-Seq data to sequences with multiple repeats and the risk of introducing false-positive variant calls, we excluded variants from the repetitive sequences in the hypervariable region 2 (HV2) of the mitochondrial control region.We determined the consensus sequence of each sample using FastaAlternateReferenceMaker (McKenna et al., 2010).
We aligned consensus mitogenome sequences from individuals with both RNA-Seq and WGS data against each other using DECIPHER R packages (Wright, 2016(Wright, , 2020) ) using default settings.
We further visualised and analysed the alignments with AliView (Larsson, 2014).

| Nuclear mitochondrial DNA segment (NUMT) identification
The presence of NUMTs in the cane toad's genome has not been documented previously.Because NUMTs can inflate estimates of genetic variations (Song et al., 2008), it is prudent to investigate their presence in mitogenome population genetic studies.We downloaded the cane toad draft genome from NCBI BioProject PRJEB24695.We used RepeatModeler (http:// www.repea tmask er.org/ Repea tMode ler/ ) to construct species-specific repeat library and masked the genome using RepeatMasker (Tarailo-Graovac & Chen, 2009).We used NUMTFinder v0.5.3 (Edwards et al., 2021) to identify potential NUMTs by searching the full and repeat-masked versions of the nuclear genome with our new mtDNA assembly.

| Population genetic analyses
We examined the spatial distribution of haplotypes by creating a median-joining haplotype network using Network and postprocessed using the maximum parsimony calculation to remove unnecessary median vectors and links (Polzin & Vahdati Daneshmand, 2003).We drew the final network using Network Publisher v2.1.2.5 (Fluxus Engineering, Clare, UK).This same procedure also was used to analyse the ND3 dataset (Table S1; Slade & Moritz, 1998).
We then investigated population genetic diversity and differentiation.We calculated the following diversity indices using DnaSP v6 (Rozas et al., 2017): nucleotide diversity (π), number of polymorphic sites, number of parsimony-informative sites, number of haplotypes (H), and haplotype diversity (h).Due to the uneven sample size across populations, we calculated haplotype richness (H R ) using FSTAT v2.9.4 (Goudet, 1995) to measure the number of alleles independent of the sample size.We used DnaSP to generate a haplotype list for Arlequin software v3.5.2.2 (Excoffier & Lischer, 2010) and a Roehl data file for use with Network software 10.2.0.0 (Bandelt et al., 1995(Bandelt et al., , 1999).To study the population structure, we employed Arlequin to calculate genetic differentiation among and between populations using pairwise fixation indices (F st ) and analysis of molecular variance (AMOVA) to estimate the genetic variation among and within populations.
To study the demographic history of each population, we used both neutrality statistics and mismatch distribution tests to infer population expansion.We produced Fu's Fs statistics implemented in Arlequin and estimated Ramos-Onsins and Rozas R 2 statistics using DnaSP.Fu's Fs test calculates the probability of observing a number of haplotypes similar to or smaller than the observed number of samples (Fu, 1997).Both tests are based on the infinite-site model without recombination, which is suitable for analysing mitochondrial genomic data.A significant negative value of Fs indicates an excess of rare alleles and haplotypes due to a sudden population growth after a population bottleneck or genetic hitchhiking, while positive values indicate a lack of rare mutations, resulting from balancing selection or population stability.Ramos-Onsins and Rozas R 2 statistic compares the difference between the number of singleton mutations and the average number of segregating sites (Ramos-Onsins & Rozas, 2002).The R 2 statistic is known to have better power than Fu's Fs when the sample size is small (<10) or the number of segregating sites is low (Ramirez-Soriano et al., 2008).A small positive value of R 2 indicates population expansion.We calculated the statistical significance of these tests using 10,000 simulations.
We constructed mismatch distributions of the frequency of pairwise differences between haplotype pairs in each population using the sudden expansion and spatial expansion models implemented in Arlequin (Excoffier, 2004;Schneider & Excoffier, 1999).We used the sum of squared deviations (SSD) to identify differences between observed and expected values and used the raggedness index to explore the smoothness of the observed mismatch distribution.
We calculated the statistical significance of these tests using 500 simulated replicates in Arlequin.A non-significant SSD value and a small raggedness index value with a unimodal distribution indicate no deviation of the observed data from the expectation under the expansion model.

| Summary of the cane toad reference mitochondrial genome
The mitogenome of R. marina (Figure 2) consists of 18,152 bp, including 22 tRNA genes, 2 rRNA genes, and 13 protein-coding genes, with 28 genes located on the forward (+) strand and 9 located on the reverse (−) strand (Table S2).There is also a 2753 bp GC-rich control region containing the origin of replication.Two groups of tandem repeats were observed in the latter part of the control region: 530 bps with 6-full repeats and 350 bps with 3-full repeats.The gene order of the cane toad mitogenome was examined by comparing it to the closely related Rhinella cf.acrolopha, and we found no gene rearrangements (data not shown).The base composition of the mitogenome was A: 29%, T: 33%, G: 14%, and C: 24%, with a bias with 62% of A + T. The codon usages of those with nucleotide G in the third codon position were lower than those with A and T. GCG (Ala) was the least used codon, and TGG was the only codon used for encoding tryptophan.All protein-coding genes used ATG as the start codon.

| Verification of mtDNA genomes from RNA-Seq data
A near-complete mitochondrial genome with 17,261 to 17,270 bps (95.1% of reference mitochondrial genome) for a total of 132 individuals was generated from the RNA-Seq and WGS dataset.A sequence with length of 889 bases consisting of repeats in the HV2 section of the control region was excluded due to difficulty in correctly mapping sequencing reads on repeat regions.To determine the accuracy of our RNA-Seq-derived mtDNA sequences, we then compared the sequences of the eight individuals for which we had both RNA-Seq and WGS data (Table S3).We identified four sites (site 2244, site 5162, site 7162, and site 7920) that showed intraindividual variation across the two datasets in multiple individuals (i.e.potential RNA editing).Two sites in coding regions (site 267 and site 477) had intra-individual differences in only one individual.Additionally, one individual had 12 sites in the control region that differed across datasets (Table S3).We replaced the sites with intra-individual variation in multiple individuals with Ns, minimising variation for subsequent population analyses that might be introduced by RNA editing.

| NUMT analysis
We found 42 short mtDNA fragments (38-140 bps in length) that mapped to 42 individual contigs of the cane toad draft nuclear genome with over 80% identity (Table S4).All fragments corresponded to the control region of the mitogenome.After the draft genome was masked for repeat elements, no mtDNA fragments were found to map to the nuclear genome.

| Population genetic analyses
We identified 262 (1.30%) polymorphic sites in the RNA-Seq dataset, 198 of which were parsimony-informative, including 25 from RNA genes, 145 from protein-coding genes, and 28 from noncoding regions and the control region (Table S2), which resulted in the identification of 58 haplotypes (Table 1; Figure 3).The 13 native F I G U R E 2 Mitochondrial genome structure of Rhinella marina.Genes are coloured according to the functional categories shown in the legend.Genes on the positive strand are on the outside of the outer circle, and genes on the negative strand are on the inside of the outer circle.The inner grey graph shows the GC content of the mitogenome.The starting and stop positions of each gene are described in Table S2.
samples contained 12 unique haplotypes.Two haplotypes found in all introduced populations accounted for nearly 50% of the sampled individuals.Within introduced populations, a "complex-star" shape topology with mostly singleton mutants was observed, and 41 of the 46 introduced haplotypes (89%) were private to a single introduction (Hawai'i: 10; Australia: 31).
We also compared the nearly complete mitochondrial genome dataset with the ND3 dataset to demonstrate the consistency of haplotype inference across studies and to put our data into the context of previous work.A total of 24 haplotypes were identified in this analysis (Figure 4): six haplotypes were identified in French Guiana/ Guyana, one haplotype was identified in Indonesian samples, three in Hawai'ian and Australian samples, and the remaining haplotypes were from other South American localities.Haplotypes from either side of the Andes clearly formed distinct lineages, which is consistent with previous work (Acevedo et al., 2016).The single Australian haplotype (Hap1) previously identified in (Slade & Moritz, 1998) was the most common haplotype in our Australian samples, and the other two haplotypes (Hap2 and Hap3) within Australia were a single nucleotide different to this common haplotype.
Haplotype diversity was highest in the native range, with reduced diversity in both Hawai'ian and Australian populations (h: 0.9782 vs. 0.8800 and 0.8344, respectively) (Table 1, Figure 5a).
Queensland toads had similar haplotype diversity but higher nucleotide diversity and haplotype richness to Hawai'ian toads.The population from the Northern Territory had noticeably lower genetic diversity and haplotype richness than other Australian populations TA B L E 1 Population genetic statistics of native-range population (French Guiana) and introduced populations (Hawai'i and Australia).

F I G U R E 3
The haplotype network illustrating the relationships among 132 cane toads (Rhinella marina) from French Guinea, Hawai'i, and Australia based on near-complete mitochondrial genomes.Each sample location is represented by a different colour.A black solid circle represents a putative haplotype that was not directly sampled.The red number adjacent to a line indicates the number of nucleotide substitutions between haplotypes, while the absence of a number indicates a single-nucleotide substitution.
(Table 1, Figure 5a,c).When the relationship between haplotypes across the invasion trajectory is considered, the 12 haplotypes from French Guiana were highly dissimilar to those from the Australian and Hawai'ian populations, forming two distinct clades (Figure 3). of 38 base pair changes (Figure 3).Australian and Hawai'ian haplotypes were intermixed, and most of the haplotypes were one base pair different to other haplotypes, with invasive haplotypes being separated from at least one other haplotype by no more than three base pair changes.
AMOVA results showed that the overall genetic variation among populations (64.49%) was larger than the variation within populations (36.33%) (Table 2).Pairwise F st values were significant between the native population and each introduced population (F st : 0.55-0.67;Mismatch distribution graphs of all populations are shown in Figure S1.The native population showed a multimodal distribution with non-significant SSD and raggedness values, indicating support for spatial and demographic expansion (Table 4).Most introduced populations showed a unimodal distribution for both expansion models with non-significant SSD and raggedness indices (Table 4), supporting the presence of spatial and demographic expansions.

| DISCUSS ION
Here, we constructed the first reference mitochondrial genome of R. marina and used nearly complete mitochondrial genome sequences to study the population genetics of the iconic cane toad in its native range and invasive ranges in Hawai'i and Australia.
We utilised RNA-Seq-derived mitogenomes to illustrate a greater number and diversity of haplotypes in mtDNA present in Australia and Hawai'i, in contrast to the surprisingly low genetic diversity previously described using the ND3 gene alone.We also found evidence of demographic and spatial expansion in invasive cane toad populations.This study provides new insights into the diversity of native cane toads and how they have evolved along their invasion trajectory.

| Considerations for data validity
The use of transcriptomic data to build mitochondrial genomes has been suggested as a cost-effective alternative to long PCRs or direct TA B L E 2 Analysis of molecular variance (AMOVA) of all populations.

Fixation indices
All populations sequencing of total DNA (Forni et al., 2019), but it is important to ensure that the inferred polymorphisms accurately reflect genetic variation.In our comparison of mtDNA data obtained from WGS and RNA-Seq, we found discrepancies in multiple nucleotide positions, all of which were in rRNA, tRNA, and the control region.These could be the result of post-transcriptional modifications or RNA editing, which are well-recognised phenomena across multiple taxa (Levanon et al., 2004;Porath et al., 2017;St Laurent et al., 2013).While common substitutions such as adenosine-to-inosine (A-to-I) deamination and cytidine-to-uridine (C-to-U) pyrimidine exchange were found, we also observed less common substitutions (e.g.thymine-to-cytidine (Tto-C)) that have been previously documented in another amphibian, Xenopus tropicalis (Zaranek et al., 2010).These processes are important for creating multiple transcript isoforms gene regulation (Tang et al., 2012) and are to be associated with adaptation (Duan et al., 2017(Duan et al., , 2021)).RNA editing of the tRNA anticodon recognition site would extend the codon-matching ability and further increase the variation of the resulting translated proteins.Although we did not observe nucleotide differences in position 34 or position 37 of tRNA, which are typically recognised for deamination of adenosine to inosine, which can be paired with A, C, or U residues (Jackman & Alfonzo, 2013), we did find changes in the D arm of tRNA (aspartic acid and lysine), the T arm of tRNA (phenylalanine), and the acceptor stem of tRNA (tryptophan) (Figure S2), which are critical for the recognition of aminoacyl tRNA synthetases to bind appropriate amino acids to its designated tRNA (Ganesh & Maerkl, 2022).Failure to form a tRNA-amino acid complex could inhibit the use of a specific amino acid in subsequent protein translation.Further proteome analysis is required to determine whether these changes affect the protein translation machinery.
Post-transcriptional modifications were more diverse in the native populations than in the invasive populations.In particular, one A search for NUMTs in the published cane toad draft genome (Edwards et al., 2018) only yielded only matches with simple repeat regions within the mtDNA control region.In other species, NUMTs originate across the entire mtDNA genome (Dayama et al., 2020;Wei et al., 2022).Following the repeat masking of the draft nuclear genome, we found no evidence of NUMTs suggesting that the putative NUMTs were just low-complexity regions in the nuclear genome.While ancient NUMTs may be too divergent to be detected by the methods employed, NUMTs with such low identity would not interfere with our mtDNA analysis.
It is worth noting that the consensus sequences generated from the bioinformatic pipelines used in this study represented the most common nucleotide at positions where multiple nucleotides were found.Although these discrepancies could result from sequencing errors, our filtering is likely to have prevented them.As discussed above, RNA editing or NUMTs could have resulted in putative polymorphisms within individuals in our RNA-Seq dataset.It is also possible that in some cases, heteroplasmy could have caused these putative polymorphisms.While determining the cause of this phenomenon is beyond the scope of this study, it is unlikely that this has affected population-level statistics.

| Population genetic analyses
Interestingly, both the haplotype network (Figure 3) and the multimodal mismatch distribution (Figure S1) indicate the presence of TA B L E 4 Demographic estimators of native-range population (French Guiana) and introduced populations (Hawai'i and Australia).multiple genetically divergent lineages within native populations.
Notably, individuals within each lineage are intermixed, rather than being grouped by geographical locations.For instance, despite an 80 km geographical separation, haplotypes H4 and H7 are found in the same lineage, while haplotypes H12 and H11 were collected only 9 m apart but belong to different lineages.Nativerange cane toads can travel 20-180 m daily for sheltering and foraging purposes (Shine et al., 2021).However, the overall dispersal rates are quite low, as the toads ultimately tend to return to the same shelters (DeVore et al., 2021).Consequently, the possibility of introgression between toads from these distinct lineages exists.
However, a recent study in South and Central America which uti-  4).In comparison with our analysis using nearly complete mitogenomes discussed above, this analysis highlights that a single gene is not sufficient to accurately infer genetic diversity.
The significant decrease in haplotype diversity and nucleotide diversity observed in both Hawai'ian and Australian populations, compared to native populations, does support that genetic bottlenecks and founder effects occurred post-introduction (Table 1, Figure 5).
This finding is consistent with previous studies on this species using nuclear genetic markers including SNPs, microsatellite markers, and MHC class I and II loci (Lillie et al., 2014(Lillie et al., , 2016;;Selechnik, Richardson, Shine, DeVore, et al., 2019).In this study, all three neutrality indices showed strong signals of population expansion and the observed data fitted to both demographic expansion and spatial expansion models in most invasive sampling sites (Table 4).This conclusion is supported by the star-shaped topology of the haplotype network created using the nearly complete mitochondrial genome dataset (Figure 3).Despite reduced diversity in the introduced range, the degree of the decrease is less than previously reported (Leblois et al., 2000).Additionally, the star-shaped topology of the network containing the introduced haplotypes, the absence of native samples here, and the significant negative Fu's Fs values indicate that many of the singleton haplotypes discovered in Hawai'i and Australia may be the products of post-introduction mutations.Alternatively, due to the limited sample sizes in the native range, it is possible that any introduced haplotypes exist but were not included in our sample set for analysis.
We observed a curvilinear pattern of genetic diversity indices across Australian populations (from Queensland to Western Australia), a pattern also found in spleen gene expression (Selechnik, Richardson, Shine, Brown, & Rollins, 2019), limb length (Hudson, McCurry, et al., 2016), spleen, and body mass (Brown, Kelehear, Shilton, Phillips, & Shine, 2015).The mechanism behind this reduction of genetic diversity in Northern Territory toads remains unclear, but it may be due to demographic processes that reduce the effective population size and increase genetic drift (Nei & Tajima, 1981), genetic hitchhiking that removes neutral polymorphic sites under positive selection (Fay & Wu, 2000), or directional selection due to harsh climate (Selechnik, Richardson, Shine, DeVore, et al., 2019), possibly in combination with dispersion-driven spatial sorting at the invasion front (Clarke et al., 2019).Despite our finding that the Northern Territory and Western Australian populations were not genetically differentiated, it is likely that they have experienced different local evolutionary forces and demographic processes.
Northern Territory toads, unlike other populations, only exhibited demographic expansion and not spatial expansion, as supported by lower levels in all statistical measures.
Our results demonstrate significant differentiation between native and introduced populations (Table 3).We identified three genetic clusters: (1) native toads, (2) Hawai'ian and Queensland toads (F st between locations: 0.02, p = .06),and (3) Northern Territory and Western Australian toads (F st between locations: 0, p = .72),consistent with previous results using SNPs from the nuclear genome (Selechnik, Richardson, Shine, DeVore, et al., 2019).It is interesting that the Queensland toads are not differentiated from those sampled in Hawai'i but differ from toads sampled in the Northern Territory and Western Australia.Queensland toads have been separated from their Hawai'ian conspecifics for almost 90 years but do share similar environmental conditions (Selechnik, Richardson, Shine, DeVore, et al., 2019).Northern Territory and Western Australian toads suffer from reduced genetic diversity compared to Queensland toads, so they have likely experienced more genetic drift, which may explain this difference.However, it is also possible that the striking reduction in precipitation and increase in temperature experienced by Northern Territory and Western Australian toads compared to those from Queensland have resulted in selective differences in the mitogenome between these populations.In support of this possibility, significant negative Fu's Fs can indicate the presence of selection.

| Conclusion
Our study used the whole mtDNA genome as a genetic marker to investigate the evolution of cane toads across their invasion trajectory from the native range to Australia, adding to our knowledge of changes to genetic diversity as cane toads have been serially introduced to multiple locations.Unsurprisingly, we found reduced genetic diversity in introduced populations, suggesting that genetic bottlenecks have occurred in these populations.Despite this reduction, our data indicate greater genetic diversity in introduced populations than has been previously described using mtDNA.We identified signatures of demographic and spatial expansion in most introduced localities we sampled, consistent with ongoing invasion.We described 46 haplotypes in introduced populations, mostly consisting of single base pair changes, suggesting evolution following introduction.By comparing WGS-and RNA-Seq-derived mitogenomes, we have identified the presence of post-transcriptional modification or RNA editing, which may supply additional genetic diversity upon which selection act.

F
Sample collection sites of the cane toad (Rhinella marina) in (a) French Guiana (the native range), (b) Hawai'i (invasive range), and (c) Australia (invasive range) for RNA-Seq and DNA-Seq.Red dots indicate samples for RNA-Seq, green squares indicate samples for DNA-Seq, and blue triangles indicate samples with both RNA-Seq and DNA-Seq.
These were separated from each other by a minimum of 110 base pair changes and from the common invasive haplotype by a minimum F I G U R E 4 Haplotype network of cane toads (Rhinella marina) from native and invasive ranges.Network built using a 510 bp segment of ND3 with flanking tRNAs.Included are a total of 223 samples, incorporating samples from this study, Slade and Moritz (1998), and Acevedo et al. (2016).Sample locations for each study are represented by different colours.A black solid circle represents a putative haplotype that was not directly sampled.The red number adjacent to a line indicates the number of nucleotide substitutions between haplotypes, while the absence of a number indicates a single-nucleotide substitution.Native-range samples are separated by the Andes Mountains in South America (Acevedo et al., 2016).Toads on the west side of the Andes are now known as Rhinella horribilis.Invasive samples are most similar to cane toads from the east of the Andes.F I G U R E 5 Haplotype diversity (a) and nucleotide diversities (b) of native population and introduced populations.
However, the Northern Territory population had significant or marginally significant p-values for SSD and raggedness values for both expansion models indicating poor support (spatial expansion: SSD pvalue = .03,raggedness p-value = .05;demographic expansion: SSD p-value = .07,raggedness p-value = .03).
individual from French Guiana displayed discrepancies in the control region between WGS and RNA-Seq data, all of which were transitions.This observation might result from tissue-specific expression in spleen, heteroplasmy, or RNA editing.Studies of RNA editing in mtDNA have primarily focused on protein-coding genes and tRNAs, with less research conducted on the control region of the mtDNA, which is responsible for regulating RNA and DNA synthesis.While the implications of potential RNA editing in cane toads remain unclear, these events may contribute to diversity at a transcriptional and translational level.NUMTs can potentially elevate estimates of mitochondrial genetic variation and contribute to false-positive signals of variation in nuclear data(Song et al., 2008).Despite NUMTs generally being perceived as non-functional due to their ability to shift reading frames, the presence of in-frame stop codons, and the use of different genetic codes, the transcription of NUMTs into non-coding RNA molecules can be captured by RNA-Seq and accidentally incorporated into assembled mitogenomes due to sequence similarity.The base differences in NUMTs will be potentially falsely interpreted as mitochondrial variation among or within individuals.Investigation of the presence/absence of NUMTs in cane toad nuclear genome and the potential expression of NUMTs is required to distinguish between genuine mtDNA variation within individuals and false-positive signals caused by incorporating nuclear-derived sequencing reads into mitochondrial population genetic analyses.
lised ddRAD sequencing has offered a different perspective on the phylogenetic relationship of cane toads.This study revealed toads from the coastal area formed two sister clades, distinct from the clade formed by toads from the rainforest regions(Mittan-Moreau et al., 2022).It becomes crucial to ascertain the precise populations of cane toads from the native range that were used as founding individuals, especially given the limited comprehensive description in prior research(Easteal, 1981).InShine et al. (2021), cane toads from the range-edge regions (Northern Territory and Western Australia) exhibited similar dispersal patterns to native toads in one of the coastal sites.This suggests the adaptation of cane toads in hot and arid environments is likely facilitated by the standing genetic variation inherited from the progenitors utilising dry, saline coastal beach environments (DeVore et al., 2021), rather than relying solely on the genetic adaptation and selection of novel mutation(s).A more thorough investigation into these distinct lineages becomes imperative for unravelling the mechanisms underlying the adaptations of cane toads in Australia.Our study challenges the assumption of low genetic diversity in Hawai'ian and Australian populations.Slade and Moritz (1998) reported a unique haplotype at the ND3 gene within Australian populations, in comparison with multiple haplotypes in native populations in Central and South America.This result supported the presence of founder effects following the introduction of cane toads into Hawai'i and Australia.With broader sampling and deeper sequencing in this study, we revealed two additional, previously undocumented, ND3 haplotypes (Hap2 and Hap3) in Hawai'ian and Australian populations, each of which differed by only one nucleotide from the previously identified haplotype (Hap1) (Figure Table3).The F st values among introduced populations were Fs values were negative in all introduced populations, indicating an excess of rare haplotypes compared to the expectation from neutrality, and all these were statistically significant.R 2 statistics were significant in Hawai'i, Queensland, and Western Australia's populations.
st : 0.01, p = .06).Within Australia, the Queensland population showed significant differentiation to both the Northern Territory (F st : 0.04, p = .01)andWesternAustralian(F st : 0.03, p = .003)populations,whereastheNorthern Territory and Western Australian populations were not genetically differentiated (F st : 0, p = .72).Two neutrality tests were implemented to study the demographic history of cane toad populations: Fu's Fs and Ramos-Onsins and Rozas' R 2 statistic (Table4).The native population had nonsignificant positive Fu's Fs values and a large R 2 statistic, indicating the population is not expanding or shrinking (i.e. is demographically stable).Fu's Population differentiation (pairwise F st ) between native range and introduced ranges (Hawai'i and Australia) (below the diagonal line).
TA B L E 3Note: The associated p-values are shown above the diagonal line.